changes: removing images0-7 and visual from cols_to_cluster¶
In [ ]:
In [ ]:
first: clone: https://github.com/AllenInstitute/visual_behavior_glm
Add to python path
In [3]:
import visual_behavior_glm
In [4]:
import visual_behavior_glm.src.GLM_params as glm_params
import visual_behavior_glm.src.GLM_analysis_tools as gat
from visual_behavior_glm.src.glm import GLM
import matplotlib.pyplot as plt
import visual_behavior.data_access.loading as loading
import visual_behavior.database as db
import plotly.express as px
import pandas as pd
import numpy as np
import os
from sklearn.cluster import KMeans
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
import seaborn as sns
In [30]:
%matplotlib inline
%widescreen
In [ ]:
get results for glm version 4¶
In [6]:
results_all = gat.retrieve_results(results_type='full')
In [7]:
results_all['glm_version'].unique()
#use v4
Out[7]:
In [8]:
rs = gat.retrieve_results(search_dict = {'glm_version': '4_L2_optimize_by_cell'}, results_type='summary')
len(rs)
Out[8]:
make a new column called 'identifier' that is {OPHYS_EXPERIMENTID}{CELL_SPECIMEN_ID}. This will give each a unique ID¶
In [9]:
rs['identifier'] = rs['ophys_experiment_id'].astype(str) + '_' + rs['cell_specimen_id'].astype(str)
In [10]:
rs
Out[10]:
pivot the table to get it into a form for clustering¶
In [11]:
rsp = rs.pivot(index='identifier',columns='dropout',values='fraction_change_from_full').reset_index()
rsp
# rsp = rs.pivot(index=['ophys_experiment_id','cell_specimen_id'],columns='dropout',values='fraction_change_from_full')
Out[11]:
merge in the identifier, cre_line, session_type and equipment_name columns¶
In [12]:
rspm = rsp.merge(rs[['identifier','cre_line','session_type','equipment_name']].drop_duplicates(),left_on='identifier',right_on='identifier',how='inner')
rspm
Out[12]:
add a column for session ID (1, 3, 4, or 6)¶
In [13]:
def map_session_types(session_type):
session_id = session_type[6:7]
return session_id
rspm['session_id'] = rspm['session_type'].map(lambda st:map_session_types(st))
rspm['session_id'].unique()
Out[13]:
save feature matrix to disk¶
In [ ]:
save = False
if save:
rspm.to_csv('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/fraction_change_var_explained_v_4_L2_fixed_lambda=1_2020.08.09.csv', index=False)
get the column names to use for clustering¶
In [18]:
cols_for_clustering = [col for col in rspm.columns if col not in ['identifier','cre_line','session_type','equipment_name', 'session_id']]
cols_for_clustering = [col for col in cols_for_clustering if col not in ['image0','image1','image2','image3','image4','image5','image6','image7',
'visual']]
cols_for_clustering
Out[18]:
do k-means on the full feature matrix¶
In [16]:
# create kmeans object
kmeans = KMeans(n_clusters=4)
In [19]:
# fit kmeans object to data
kmeans.fit(rspm[cols_for_clustering])
Out[19]:
In [20]:
rspm['k-means_cluster_id'] = kmeans.fit_predict(rspm[cols_for_clustering])
it is assigning almost all cells to a single cluster - not useful!¶
In [21]:
rspm['k-means_cluster_id'].value_counts()
Out[21]:
make a tsne plot¶
In [22]:
tsne = TSNE(n_components=2, verbose=1, perplexity=40, n_iter=300)
In [23]:
tsne_results = tsne.fit_transform(rspm[cols_for_clustering])
assign tsne reults to new columns¶
In [24]:
rspm['tsne-2d-one'] = tsne_results[:,0]
rspm['tsne-2d-two'] = tsne_results[:,1]
In [25]:
rspm['session_type'].unique()
Out[25]:
make a histogram of the omission dropout values¶
In [26]:
fig,ax=plt.subplots()
ax.hist(rspm['omissions'],bins=np.arange(-1,1,0.01))
ax.axvline(0,color='black')
ax.set_xlabel('omissions dropout fractional change in var explained')
Out[26]:
make a new column that identifies cells as 'modulated' by running or omissions (the threshold is up for debate!)¶
In [27]:
def is_modulated(v, threshold=0):
return v < threshold
rspm['omission_modulated'] = rspm['omissions'].map(lambda v:is_modulated(v, threshold=0))
rspm['running_modulated'] = rspm['running'].map(lambda v:is_modulated(v, threshold=0))
plot the tsne results, coloring dots in different ways¶
It looks like they cluster almost perfectly by omission modulation!
In [32]:
fig,ax=plt.subplots(2,2,figsize=(15,12),sharey=True,sharex=True)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="cre_line",
palette=sns.color_palette("hls", 3),
data=rspm,
legend="full",
alpha=0.3,
ax=ax[0,0]
)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_type",
palette=sns.color_palette("hls", 8),
data=rspm,
legend="full",
alpha=0.3,
ax=ax[0,1]
)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="k-means_cluster_id",
palette=sns.color_palette("hls", 4),
data=rspm,
legend="full",
alpha=0.3,
ax=ax[1,0]
)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue='omission_modulated',
palette=sns.color_palette("hls", 2),
data=rspm,
legend="full",
alpha=0.3,
ax=ax[1,1]
)
fig.tight_layout()
for ii,hue in enumerate(['cre_line','session_type','k-means label', 'omission modulated']):
ax.flatten()[ii].set_title('color by {}'.format(hue))
fig.tight_layout()
plot tsne again, grouped by cre-line¶
In [33]:
fig,ax=plt.subplots(2,2,figsize=(20,12),sharey=True,sharex=True)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm, #.query('cre_line=="Vip-IRES-Cre"'),
legend='full',
alpha=0.3,
ax=ax[0,0]
)
ax[0,0].set_title('all cre lines')
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm.query('cre_line=="Vip-IRES-Cre"'),
legend='full',
alpha=0.3,
ax=ax[0,1]
)
ax[0,1].set_title('VIP only')
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm.query('cre_line=="Slc17a7-IRES2-Cre"'),
legend='full',
alpha=0.3,
ax=ax[1,0]
)
ax[1,0].set_title('Slc17a7 only')
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm.query('cre_line=="Sst-IRES-Cre"'),
legend='full',
alpha=0.3,
ax=ax[1,1]
)
ax[1,1].set_title('Sst only')
Out[33]:
plot tsne again, colored continuously by the omissions dropout¶
In [34]:
fig,ax=plt.subplots(1,1,figsize=(10,6),sharey=True,sharex=True)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue='omissions',
palette='magma',
data=rspm.query('omissions>-1 and omissions<1'),
legend=False,
alpha=0.3,
ax=ax
)
Out[34]:
do PCA¶
a huge amount of variance is in just the first few PCs!
In [35]:
rspm
Out[35]:
In [36]:
cols_for_clustering
Out[36]:
In [37]:
np.shape(rspm[cols_for_clustering].values)
Out[37]:
In [38]:
pca = PCA(n_components=len(cols_for_clustering))
pca_result = pca.fit_transform(rspm[cols_for_clustering].values)
rspm['pc1'] = pca_result[:,0]
rspm['pc2'] = pca_result[:,1]
rspm['pc3'] = pca_result[:,2]
print('Explained variation per principal component: {}'.format(pca.explained_variance_ratio_))
In [40]:
np.shape(pca.components_)
Out[40]:
In [41]:
fig,ax=plt.subplots()
ax.plot(
np.arange(len(cols_for_clustering)),
pca.explained_variance_ratio_,
'o-k'
)
ax.set_xlabel('PC number')
ax.set_ylabel('variance explained')
Out[41]:
In [42]:
np.shape(pca_result)
Out[42]:
In [43]:
np.shape(pca_result[0,:])
Out[43]:
examine the PC components individually¶
followups
- remove redundant features before performing PCA
- Only plot PC dimensions that explain a reasonable amount of variance (each PC will sum to one, so low variance PCs will be misleading).
In [44]:
pca.components_[:,0].sum()
Out[44]:
In [45]:
pca.components_[:,0]
Out[45]:
In [46]:
fig,ax=plt.subplots(figsize=(15,5))
N_PCs = 10
for PC in range(N_PCs):
ax.plot(
pca.components_[:,PC]
)
ax.legend(np.arange(N_PCs))
ax.axhspan(1/np.sqrt(32),-1/np.sqrt(32),zorder=-np.inf,alpha=0.25)
ax.set_xticks(np.arange(32))
ax.set_xticklabels(cols_for_clustering,rotation=45,ha='right')
fig.tight_layout()
Can I do k-means in PCA space?¶
No, it is still assigning nearly all cells to a single cluster
In [47]:
pc_cols = ['pc1','pc2','pc3']
kmeans = KMeans(n_clusters=3)
kmeans.fit(rspm[pc_cols])
rspm['k-means_on_PCs'] = kmeans.fit_predict(rspm[pc_cols])
rspm['k-means_on_PCs'].value_counts()
Out[47]:
plot PCA results, colored by cre-line¶
In [48]:
fig,ax = plt.subplots(figsize=(16,10))
sns.scatterplot(
x="pc1", y="pc2",
hue="cre_line",
palette=sns.color_palette("hls", 3),
data=rspm,
legend="full",
alpha=0.3,
ax=ax
)
ax.set_xlim(-1,1)
ax.set_ylim(-1,1)
Out[48]:
In [ ]:
rspm.query('cre_line=="Sst-IRES-Cre"')
In [70]:
rspm.keys()
Out[70]:
In [53]:
sample_size=1000
dfs = []
for cre_line in rspm['cre_line'].unique():
dfs.append(rspm.query('cre_line == @cre_line').sample(sample_size, random_state=0))
data_sample = pd.concat(dfs)
data_sample
Out[53]:
In [54]:
data_sample.keys()
Out[54]:
make a 3D plotly plot of the PCA results (1000 cells from each cre line)¶
In [52]:
query_string = '''pc1>-1 and pc1 < 1 and pc2>-1 and pc2 < 1 and pc3>-1 and pc3 < 1'''
fig = px.scatter_3d(
data_sample.query(query_string),
x='pc1',
y='pc2',
z='pc3',
color='cre_line',
)
fig.update_traces(
marker=dict(
size=3,
opacity=0.25
)
)
fig.update_layout(
margin=dict(l=30, r=30, t=10, b=10),
width=1200,
height=1000,
)
# fig.write_html("/home/dougo/code/dougollerenshaw.github.io/figures_to_share/2020.08.09_PCA_on_GLM.html")
# fig.show()
do kmeans clustring on the tsne space¶
It kinda works
In [57]:
kmeans = KMeans(n_clusters=2)
kmeans.fit(rspm[pc_cols])
tsne_cols = ['tsne-2d-one','tsne-2d-two']
rspm['k-means_on_tSNE'] = kmeans.fit_predict(rspm[tsne_cols])
rspm['k-means_on_tSNE'].value_counts()
Out[57]:
In [58]:
fig,ax=plt.subplots(1,1,figsize=(10,6),sharey=True,sharex=True)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue='k-means_on_tSNE',
palette=sns.color_palette("hls", 2),
data=rspm,
legend='full',
alpha=0.3,
ax=ax,
)
Out[58]:
try agglomerative clustering?¶
It's also assigning to a single cluster
In [59]:
# import hierarchical clustering libraries
import scipy.cluster.hierarchy as sch
from sklearn.cluster import AgglomerativeClustering
In [60]:
# create clusters
hc = AgglomerativeClustering(n_clusters=4, affinity = 'euclidean', linkage = 'ward')
# save clusters for chart
y_hc = hc.fit_predict(rspm[cols_for_clustering])
In [61]:
rspm['hc'] = y_hc
rspm['hc'].value_counts()
Out[61]:
In [62]:
import umap
reducer_2d = umap.UMAP(n_components=2)
reducer_3d = umap.UMAP(n_components=3)
embedding_2d = reducer_2d.fit_transform(rspm[cols_for_clustering])
embedding_3d = reducer_3d.fit_transform(rspm[cols_for_clustering])
embedding_3d.shape
Out[62]:
In [ ]:
assign umap results to new columns¶
In [63]:
rspm['umap_3d_embedding_0'] = embedding_3d[:, 0]
rspm['umap_3d_embedding_1'] = embedding_3d[:, 1]
rspm['umap_3d_embedding_2'] = embedding_3d[:, 2]
rspm['umap_2d_embedding_0'] = embedding_2d[:, 0]
rspm['umap_2d_embedding_1'] = embedding_2d[:, 1]
do kmeans on the 2d embedding¶
only two clusters since the 2d umap gives two clear clusters k-means seems to identify them perfectly!
In [64]:
kmeans = KMeans(n_clusters=2)
umap_2d_cols = ['umap_2d_embedding_0','umap_2d_embedding_1']
rspm['k-means_on_2d_umap'] = kmeans.fit_predict(rspm[umap_2d_cols])
rspm['k-means_on_2d_umap'].value_counts()
Out[64]:
do kmeans on the 3d embedding¶
In [65]:
kmeans = KMeans(n_clusters=3)
umap_cols = ['umap_3d_embedding_0','umap_3d_embedding_1','umap_3d_embedding_2']
rspm['k-means_on_3d_umap'] = kmeans.fit_predict(rspm[umap_cols])
rspm['k-means_on_3d_umap'].value_counts()
Out[65]:
plot the 3d results in 2d¶
k-means clearly didn't do a good job of identifying two of the three clusters
In [66]:
fig,ax=plt.subplots(1,1,figsize=(10,6),sharey=True,sharex=True)
sns.scatterplot(
x='umap_3d_embedding_0',
y='umap_3d_embedding_1',
hue='k-means_on_3d_umap',
palette=sns.color_palette("hls", 3),
data=rspm,
legend='full',
alpha=0.3,
ax=ax,
)
Out[66]:
In [71]:
import seaborn as sns
fig,ax=plt.subplots(2,2,figsize=(15,12),sharey=True,sharex=True)
sns.scatterplot(
x='umap_3d_embedding_0',
y='umap_3d_embedding_1',
hue="cre_line",
palette=sns.color_palette("hls", 3),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[0,0]
)
sns.scatterplot(
x='umap_3d_embedding_0',
y='umap_3d_embedding_1',
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[0,1]
)
sns.scatterplot(
x='umap_3d_embedding_0',
y='umap_3d_embedding_1',
hue='k-means_on_2d_umap',
palette=sns.color_palette("hls", len(rspm['k-means_on_2d_umap'].unique())),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[1,0]
)
sns.scatterplot(
x='umap_3d_embedding_0',
y='umap_3d_embedding_1',
hue='omission_modulated',
palette=sns.color_palette("hls", 2),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[1,1]
)
fig.tight_layout()
for ii,hue in enumerate(['cre_line','session_type','k-means label', 'omission modulated']):
ax.flatten()[ii].set_title('color by {}'.format(hue))
fig.tight_layout()
now make a 3d plotly plot of the same¶
In [67]:
import plotly.express as px
fig = px.scatter_3d(
rspm,
x='umap_3d_embedding_0',
y='umap_3d_embedding_1',
z='umap_3d_embedding_2',
color='k-means_on_3d_umap'
)
fig.show()
plot the 2d umap results¶
In [68]:
import seaborn as sns
fig,ax=plt.subplots(2,2,figsize=(15,12),sharey=True,sharex=True)
sns.scatterplot(
x='umap_2d_embedding_0',
y='umap_2d_embedding_1',
hue="cre_line",
palette=sns.color_palette("hls", 3),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[0,0]
)
sns.scatterplot(
x='umap_2d_embedding_0',
y='umap_2d_embedding_1',
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[0,1]
)
sns.scatterplot(
x='umap_2d_embedding_0',
y='umap_2d_embedding_1',
hue='k-means_on_2d_umap',
palette=sns.color_palette("hls", len(rspm['k-means_on_2d_umap'].unique())),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[1,0]
)
sns.scatterplot(
x='umap_2d_embedding_0',
y='umap_2d_embedding_1',
hue='omission_modulated',
palette=sns.color_palette("hls", 2),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[1,1]
)
fig.tight_layout()
for ii,hue in enumerate(['cre_line','session_type','k-means label', 'omission modulated']):
ax.flatten()[ii].set_title('color by {}'.format(hue))
fig.tight_layout()
plot umap with omissions modulation as continuous color¶
In [69]:
fig,ax=plt.subplots(1,1,figsize=(10,6),sharey=True,sharex=True)
sns.scatterplot(
x='umap_2d_embedding_0',
y='umap_2d_embedding_1',
hue='omissions',
palette='magma',
data=rspm.query('omissions>-1 and omissions<0.25'),
legend=False,
alpha=0.3,
ax=ax,
)
Out[69]:
In [ ]: